
program main

    use solve, only: solve1, Gauss, solve2
    real(kind=8), allocatable :: A(:, :), b(:), c_(:), c1(:), b_(:, :)
    real :: tic, toc
    real :: time1, time2, time3

    integer, parameter :: N = 3000
    integer :: ipiv(N)

    allocate (A(N, N), b(N), c_(N), c1(N), b_(N, 1))

    call random_number(A)
    call random_number(b)

    b_(:, 1) = b
    ! call cpu_time(tic)
    ! call DGESV( N, 1, A, N, IPIV, b_, N, INFO )
    ! call cpu_time(toc)
    ! print "(A, ES12.4)", "solution4(LAPACK)      , toc - tic(seconds): ", toc - tic
    call cpu_time(tic)
    b_ = solve1(A, b_)
    call cpu_time(toc)
    time1 = toc - tic
    print "(A, ES12.4)", "solution6(Euler2)      , toc - tic(seconds): ", time1

    deallocate (b_)
    allocate (b_(N, N + 1))
    b_(:, 1:N) = A
    b_(:, 1) = b
    
    call cpu_time(tic)
    call solve2(b_, N, 1)
    call cpu_time(toc)
    time2 = toc - tic
    print "(A, ES12.4)", "solution6(Euler3)      , toc - tic(seconds): ", time2
    
    call cpu_time(tic)
    call Gauss(b_, N, N + 1)
    call cpu_time(toc)
    time3 = toc - tic
    print "(A, ES12.4)", "solution*(Walcs )      , toc - tic(seconds): ", time3
    
    print *, "solution4(Walcs ) / solution6(Euler) = ", time3 / time2

    deallocate (a, b, c_, c1, b_)
    allocate (a(2, 2), b(2), b_(2, 3))
    a = reshape([1, 3, 2, 4], [2, 2])
    b = [1, 1]
    b_(:, 1:2) = A
    b_(:, 3) = b
    call Gauss(b_, 2, 3)
    print *, b_(:, 3)
    print *, solve1(a, reshape(b, [2, 1]))

end program main